Towards Optimizing Explicit Fracture Simulations of Unconventional Reservoirs

Daniel Barreca

Summary

Simulations are carried out using MRST’s black oil module on reservoirs with properties similar to the Bakken Shale. Simulations of a fractured reservoir based on fracture simulations are also carried out. Hydraulic fractures are treated through the explicit fracture modelling method. In an attempt to reduce grid complexity, the effect of fracture geometry and fracture interaction is first investigated. Results demonstrated that production is dependent on fracture geometry, therefore the mesh cannot be simplified to structured meshed for explicit fracture systems. Additionally, the method of snapshots is applied to one of the reservoirs. Results from computing the basis functions and eigenvalues from this snapshot matrix demonstrate the potential of reduced order modelling on explicitly fractured unconventional simulations.

Introduction

Due to the large amount of oil and gas held in unconventional reservoirs, it is important that efficient and accurate models exist for simulation of these resources. When trying to model unconventional reservoirs, certain challenges arise such as matrix/fracture interactions, nonlinear pressure distributions, fracture networks, and geomechanical considerations (Farah 2016). Modelling complex physical processes such as these may be computationally demanding, facilitating a need for methods that reduce CPU time of simulation. One way of reducing complexity of a problem is through a reduced order model (ROM).

Reduced order modelling is used to reduce computational demands of a simulation by reducing the degrees of freedom of a problem. Various ROM techniques exist, though the focus of this project is using a classification of ROMs known as Proper Orthogonal Decomposition (POD). More specifically, to attempt to reduce the computational demands of unconventional reservoir simulation, the goal of this project is to test the applicability of the method of snapshots to such reservoirs. The method of snapshots, developed by Sirovich, is a technique of producing a POD basis for a ROM (Sirovich 1987). POD methods have been used to create ROMs for simple reservoir simulations in a number of studies (Jing et al. 2017, Cardoso et al. 2009, Markovinovic et al. 2002). Other POD studies more specifically related to unconventional simulation include ROMs for geomechanics and coupled flow (Florez and Gildin 2018), multiscale fracture simulation (Oliver et al. 2017), and thermal simulation of oil shale (Rousset 2010). For this project, explicitly fractured simulations are first simulated and compared. The method of snapshots will then be applied explicitly fractured systems, as well as be applicable to simulation of my chosen reservoir (Bakken).

Hydraulic Fracturing

Reservoirs such as shale have extremely low values of porosity and permeability. Hydraulic fracturing is a technique where fluid is injected into a rock formation and proppants are embedded, allowing for a highly permeable flow path towards the wellbore. This method has made unconventional resources an attractive source of economic value. As a result, simulation of these resources (along with hydraulic fractures), must be accurately calculated. As one of the steps in building the reservoir model, hydraulic fracture pathway must be considered since this fracture will need to be modelled.

Though many models assume the fracture will extend linearly, the hydraulic fracture geometry may not always extend perpendicular and linearly from the wellbore. Various factors can generate more complex fracture geometries. For example, in cases where the reservoir is naturally fractured, the hydraulic fracture may follow along this pathway. Another consideration for fracture geometry is the stress field within the reservoir. There is a possibility of hydraulic fractures interacting with each other due to created stresses and depending on the fracture spacing (Wu 2014). Additionally, fracture path may be influenced by the pre-existing stress field, as indicated in fracture simulations by Wu (2014), shown in Figure 1. Hydraulic fracturing possibilities such as these will be investigated through reservoir simulation.

Preprocessing Procedure for Simulating More Complex Geometries

Because of the geometric complexity of hydraulic fractures, more advanced meshing and geometry generation software packages are used for preprocessing. Autodesk INVENTOR is a commercial 3D CAD software that is commonly used in various industries for CAD creation. The software package is used here due to its flexibility and advanced geometry creation features. Various fracture models can be created, including complex fracture geometries. A simple example is shown in Figure 2.

The CAD file is exported as a (.dwg) file and imported into ANSYS ICEM CFD. ICEM CFD is a commercial mesh generation software package. Again, a different software is used here for meshing because ICEM CFD is more flexible than MATLAB in importing CAD files. In addition to this, ICEM CFD has more advanced meshing control and generation options than MRST. Cell sizings, cell growth rates, and various meshing algorithms can easily be implemented to mesh the geometry. An example of a generated triangular mesh is again shown in Figure 2. Cell sizing can be straightforwardly refined or coarsened by specifying a cell size or number of nodes around specified curves, allowing for controlled node density and growth rates around hydraulic fractures.

From here, the mesh can be exported into specific formats for various simulation packages such as FLOW3D or ANSYS FLUENT. MRST is not one of the available options for mesh exporting. As a workaround, the mesh is exported as a (.msh) text file, which is the unstructured mesh format that ANSYS FLUENT accepts. X and Y coordinates from the nodal positions are extracted from this text file and stored as vectors in MATLAB. Finally, using these two vectors, the mesh can be reconstructed using MRST’s TriangularGrid() function.

Fracture permeability and porosity cannot simply be assigned the correct values based on location because fracture geometry may be complex. Instead, to assign permeability and porosity to the fracture, all cells are looped over and the volumes are checked. A histogram can be created to determine which size cells lie within the fracture. Since the mesh is most refined in these sections, cell volumes should be smallest within the fractures. Unfortunately, this is not always the case. As shown in Figure 2, using this method will result in slight deviations from the true fracture, as random cells outside, but near, the fracture may be included. However, a decent representation of the fracture using this method is attained. A more accurate method of assigning correct permeability values would include some implementation where cell centroid locations are utilized.

Fig. 2Preprocessing Procedure

General Simulation and Geometry Details

Because the Bakken shale is an oil shale, the simulations are carried out using MRST’s black oil module. All simulations are completed using a 2D discretization due to the computational requirements of the 3D meshes (further explained in a later section). The geometries for all following simulations have the dimensions of one of the figures specified in Figure 4. For the simulations with 3 hydraulic fractures, a fracture spacing of 250 ft is used, which is based off common spacings from fracturing done in the Bakken Shale. For the simulations consisting of 4 hydraulic fractures, a spacing of 150 ft is sued, which is based off simulations from Wu (2014). Most reservoir and fluid properties are based off the Bakken Shale and specified in the literature review. Properties not clearly defined or not included in the literature review are included in Figure 3.

Property	Value	Unit
Matrix Porosity	6	%
Matrix Permeability	0.01	mD
Fracture Porosity	20	%
Fracture Permeability	10	mD
Fracture Width	0.2	in
Wellbore Radius	0.4	in
Wellbore Pressure	1000	psi
Simulation Time	10	years
Simulation Timestep	0.0025	years

Fig. 3Extra Details for Simulation, not included in Literature Review

Fig. 4 Base Geometry Dimensions

Grid Independence Study

To ensure special convergence of the simulations, a grid independence study is initially carried out. Because pressure drawdowns, fluid properties, and rock properties are similar between all geometries, only one grid independence study is completed. The grid independence study is run this single geometry, and a similar global mesh sizing is to be used for the rest of the geometries. Three meshes of different node densities are created for the straight-line perforations’ geometry. 1000, 500, and 250-feet cell-sizings are used for the global mesh size, resulting in 14000, 7000, and 5000 nodes (Fig. 5). The primary value of interest, oil produced, is monitored and recorded over a 10 year simulation period. The results are displayed in Figure 6. While a significant difference in production exists between the 5000 node and the 7000 node system, this difference is insignificant between the 7000 and 14000 node systems. As a result, a global grid spacing of 500 ft will be specified in all following simulations.

Fig. 5Grid Independence Study: 5000 nodes (left), 7000 nodes (middle), and 14000 nodes (right)

Fig. 6Grid Independence Study, Monitoring Cumulative Oil Production

Impact of Fracture Geometry

Before moving on, an initial and important question is whether fracture geometry will have a significant impact on production. If fracture geometry does not have an impact on production, then there is little point in modelling anything beyond simple, straight hydraulic fractures. Straight hydraulic fractures have an inherent advantage for simulation in that they can be modelled easily on structured meshes. Structured meshes have various advantages such as computational efficiency and alignment with flow or permeability fields. Setting up a structured mesh for a more complex hydraulic fracture system (such as that in Figure 2) would be at the least very difficult to set up and a manually labor-intensive process. In many cases, structured meshes may even be practically impossible to generate.

Two situations are tested in this study: effect of fracture curvature (Fig. 7) and effect of fracture interaction (Fig. 8). For the fractures in Figure 7, it was ensured during the geometry creation process that every fracture had a half-length of 400 ft, spacing of 250 ft, and node density throughout the mesh. Thus, all fractures have an equivalent pore volume, and the only difference is the degree of curvature. For the fractures in Figure 8, the effect of fracture interaction is simulated by roughly replicating a geometry from Wu (2014).

Simulations are run for 10 years at a timestep of 0.0025 years with a wellbore pressure of 1000 psi. Cumulative production is plotted for each situation in figures 9 and 10. Additionally, pressure contour plots are generated and displayed in figures 11 and 12. The results demonstrate that fracture geometry may have an impact on cumulative production.

In the case of analyzing fracture curvature, a significant curvature decreases production, while a slight curvature increases production. The contour plots suggest that an increase in curvature results in an increase in near wellbore production. This can be seen at a timestep of 1 year, where pressure in areas surrounding the wellbore is lower for higher curvature fractures. This suggests that there may be some optimal curvature or angle that optimizes production. Further analysis would need to be completed to fully investigate this possibility.

In the case of analyzing fracture interaction, the results are as expected. When considering fracture interaction, the system with shortened inner fractures results in a noticeably lower production value. As demonstrated in the contour plots, this is simply due to longer fractures permeating further into the reservoir, and consequently draining more of the reservoir. One factor that was not considered here was that inner fractures merging with the main fractures may extend the main fractures’ length. Further literature review needs to be completed on this topic to fully understand hydraulic fracture interaction.

Fig. 7:Varying Hydraulic Fracture Curvature

Fig. 8Varying Hydraulic Fracture Interaction

Fig. 9Cumulative Production Plot of Results from Simulations in Fig. 7

Fig. 10Cumulative Production Plot of Results from Simulations in Fig. 8

Fig. 11Pressure (psi) Contour Plots for 5 Time Snapshots: Effect of Curvature

Fig. 12Pressure (psi) Contour Plots for 5 Time Snapshots: Effect of Fracture Interaction

Reduced Order Modelling in Reservoir Simulation

While simpler reservoir and fracture simulations such as these may not require a large amount of computational resources, situations such as 3D simulations, realistic sized reservoirs, and increased number of fractures may take a significant amount of time to complete. Take Figure 13 as an example. When no fracture is modelled, the mesh only consists of 710 nodes. After a fracture is added in and discretized with the reservoir, the node count drastically increases to 44,913 nodes, which is 63 times the number of nodes in the standard simulation. This also does not take into consideration the mesh shown is a 2D system. When extending to 3D, it may not be as simply as creating a few layers. Since the cell sizes around the fractures are on the order of magnitude of inches, the vertical resolution should also be this order of magnitude in order to maintain favorable aspect ratios and other important cell quality criteria.

Fig. 13Discretization Increase When Including Fracture

In addition to a drastic increase in cell count, an engineer may want to run hundreds or thousands of simulations related to production optimization or sensitivity analysis. The question proposed is how this high dimensional system can be solved in a more efficient manner. A potential solution to this question is reduced order modelling. There are various methods of reduced order modelling, one of those being Proper Orthogonal Decomposition (POD) based methods.

The time evolution of pressure depends on a nonlinear function of variables such as space, time, and rock properties.

Normally, using MRST, a finite volume discretization of this system will be carried out, which is a finite representation of the infinite dimensional space. Essentially, the PDE is substituted for a large system of solvable differential equations. Because explicit fracture systems significantly increase the discretization of the reservoir system, we would like to reduce this large system in some way. The first step in the reduced order modelling process is to assume that the spatial (x) dependence and the time (t) dependence of a problem can be separated through the separation of variables:

The functions of interest that need to be determined in the equation above are the basis functions, . Afterwards, a finite sum of basis functions may be determined based on the time dynamics of the system.

To determine these basis functions, the full MRST simulation must be run for a number of time steps so that a snapshot matrix for the variables of interest can be attainted. For example, the snapshot matrix, P, is simply just a collection of the pressure values of each grid cell at time interval ‘NT’. This information is stored in an (n x NT) matrix where n is the number of grid cells.

After this matrix is stored, an SVD may be computed on the matrix set.

The matrix U, which is computed as part of SVD, contains information about the dominant modes for the basis functions.

The process described above is completed for the fracture interaction system from Figure 13. Eigenvalues are plotted in Figure 14 and basis function energy contour plots are shown in Figure 15. This system may be suitable for reduced order modelling, as indicated by Figure 14. The eigenvalues in figure 14 demonstrated that the energy of the system varies by orders of magnitudes. Therefore, the vast majority of energy in this system can be represented by a small number of modes. The energy of these modes can be viewed and compared, shown in Figure 15. The first 8 basis functions(which are the most dominant energy systems) are vastly different than, for example, the 200^th largest energy mode (which is essentially just noise).

Fig. 14Semilog plot of Eigenvalues for the Fracture Interaction System

Fig. 15Basis Functions of the Complex Fracture Simulation

Conclusions

1.) Modelling accurate fracture geometries during the simulation process is important for accurate estimation of resources, Fracture geometry may have an impact on production. In some cases, fracture curvature can increase production, though a large curvature may significantly decrease production. Hydraulic fracture interaction should also be considered.

2.) Modelling fractures through an explicit method will significantly increase cell count. Reduced order modelling of unconventional has potential, as indicated by an SVD on the system.

Future Work

1.) More work needs to be done in determining effect of fracture geometry. A larger number of fractures, different fracture spacings, various reservoir properties or fluid properties, and fracture lengths in comparison to reservoir length should be investigated to fully understand the impact of fracture geometry.

2.) Different scenarios of fracture interaction may be considered to estimate the impact on production

3.) Reduced order modelling (either non-intrusive or intrusive) should be carried out. The results in this study simply show the potential of reduced order modelling, but a working reduced order model was never complete due to possible errors in the code or lack of understanding of the problem.

References

Farah, N., 2016. Flow Modelling in Low Permeability Unconventional Reservoirs (Doctoral dissertation, Paris 6).

Sirovich, L., 1987. Turbulence and the dynamics of coherent structures. I. Coherent structures. Quarterly of applied mathematics, 45(3), pp.561-571.

Cao Jing, Zhao Hui, Yu Gaoming, and Xie Yunhong, 2017. The POD Model Order Reduction Method of Accelerating Reservoir Numerical Simulation. IOSR Journal of Applied Geology and Geophysics, Volume 5, Issue 1 Ver. I, pp. 107-114.

Cardoso, M.A., Durlofsky, L.J. and Sarma, P., 2009. Development and application of reduced‐order modeling procedures for subsurface flow simulation. International journal for numerical methods in engineering, 77(9), pp.1322-1350.

Markovinovic, R., Geurtsen, E.L., Heijn, T. and Jansen, J.D., 2002, September. Generation of Low-Order Reservoir Models Using POD, Emperical Grammians and Subspace Identification. In ECMOR VIII-8th European Conference on the Mathematics of Oil Recovery.

Florez, H. and Gildin, E., 2018, September. Model-Order Reduction Applied To Coupled Flow and Geomechanics. In ECMOR XVI-16th European Conference on the Mathematics of Oil Recovery.

Oliver, J., Caicedo, M., Huespe, A.E., Hernández, J.A. and Roubin, E., 2017. Reduced order modeling strategies for computational multiscale fracture. Computer Methods in Applied Mechanics and Engineering, 313, pp.560-595.

Rousset, M., 2010. Reduced-order modeling for thermal simulation (Doctoral dissertation, Stanford University).

Wu, K., 2014. Numerical modeling of complex hydraulic fracture development in unconventional reservoirs (Doctoral dissertation).